
type HybridTaus
	z1::Uint
	z2::Uint
	z3::Uint
	z4::Uint
end


# S1, S2, S3, and M are all constants, and z is part of the generator state.  
tausStep!(z::Uint, S1::Uint, S2::Uint, S3::Uint, M::Uint) = 
	(z = (((z & M) << S3) $ (((z << S1) $ z) >> S2)))

lcgStep!(z::Uint, A::Uint, C::Uint) = (z = A * z + C)

# combined period is lcm(p1,p2,p3,p4)~ 2^121  
rng(t::HybridTaus) = 
	tausStep!(t.z1, uint(13), uint(19), uint(12), uint(4294967294)) $ 
		tausStep!(t.z2, uint(2), uint(25), uint(4), uint(4294967288)) $ 
		tausStep!(t.z3, uint(3), uint(11), uint(17), uint(4294967280)) $ 
		lcgStep!(t.z4, uint(1664525), uint(1013904223)) 

function remove!{CONT, EL}(cont::CONT, el::EL)
	for i in cont
		if i == e
			i = last(cont)
			pop!(cont)
			return
		end
	end
end


rng_state = HybridTaus(uint(5093), uint(384938), uint(93010), uint(34384))
rng() = rng(rng_state)

type Individual
	x::Int
	y::Int
	count::Int
end

type Patch
	population::Vector{Individual}

	max_x::Int
	max_y::Int

	cell_size::Int
	cache::Matrix{Vector{Int}}

	visibility::Int
	max_sq_dist::Int
end

cache_pos(patch::Patch, x::Int) = fld((x-1), patch.cell_size) + 1

torus(max::Int, x::Int) = (x+max-1)%max + 1

function Patch(pop_size, max_x, max_y, visibility, cell_size)
		
	pop = Individual[]
	for i in 1:pop_size
		push!(pop, Individual(int(rng() % max_x + 1), int(rng() % max_y + 1), 0))
	end

	cache = Array(Vector{Int}, div(max_x, cell_size), div(max_y, cell_size))
	for x in 1:size(cache, 1)
		for y in 1:size(cache, 2)
			cache[x, y] = Int[]
		end
	end

	patch = Patch(pop, max_x, max_y, cell_size, cache,
		visibility, visibility * visibility)

	for i in 1:length(patch.population)
		ind = patch.population[i]
		push!(patch.cache[cache_pos(patch, ind.x), cache_pos(patch, ind.y)], i)
	end

	return patch
end

function step(patch::Patch)
	move(patch)
	count(patch)
end

function move(patch::Patch)
	i = 1
	for ind in patch.population

		const cp_x = cache_pos(patch, ind.x)
		const cp_y = cache_pos(patch, ind.y)

		const new_x::Int = torus(patch.max_x, ind.x + int(rng()%5) - 2)
		const new_y::Int = torus(patch.max_y, ind.y + int(rng()%5) - 2)

		const cpn_x = cache_pos(patch, new_x)
		const cpn_y = cache_pos(patch, new_y)

		if cpn_x != cp_x || cpn_y != cp_y
			remove!(patch.cache[cp_x, cp_y], i)
			push!(patch.cache[cpn_x, cpn_y], i)
		end	
		
		ind.x = new_x
		ind.y = new_y

		i += 1
	end
end

function count(patch::Patch)
	const cache_m_x = size(patch.cache, 1)
	const cache_m_y = size(patch.cache, 2)

	for ind in patch.population

		const min_cp_x = cache_pos(patch, ind.x - patch.visibility)
		const min_cp_y = cache_pos(patch, ind.y - patch.visibility)
		const max_cp_x = cache_pos(patch, ind.x + patch.visibility)
		const max_cp_y = cache_pos(patch, ind.y + patch.visibility)

		for cpx in min_cp_x:max_cp_x
		
			# remember whether we are beyond the boundary
			const offs_x = cpx<1 ? -patch.max_x : (cpx>cache_m_x ? patch.max_x : 0)
			# but we also need the actual grid index
			const rx = torus(cache_m_x, cpx)

			for cpy in min_cp_y:max_cp_y
				
				const offs_y::Int = 
					cpy<1 ? -patch.max_y : (cpy>cache_m_y ? patch.max_y : 0)
				const ry = torus(cache_m_y, cpy)

				@inbounds const loc_p = patch.cache[rx, ry]
				for l::Int in loc_p
					# measure distance to projected position 
					# (needed for areas around edges)
					@inbounds const dx::Int = ind.x-(patch.population[l].x + offs_x)
					@inbounds const dy::Int = ind.y-(patch.population[l].y + offs_y)
					const sq_dist = dx*dx + dy*dy
					if sq_dist < patch.max_sq_dist
						ind.count+=1
					end
				end
			end
		end
	end
end


function run(steps::Int)

	world = [ Patch(256, 1024, 1024, 10, 64) for i in 1:256 ]

	for s in 1:steps
		println(s)
		for p in world
			step(p)
		end
	end
end	
